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Abstract 

Supernova remnants are likely to be the accelerators of the galactic cosmic rays. Assuming the 
correctness of this hypothesis, we develop a method to extract the parent cosmic ray spectrum from 
the VHE gamma ray flux emitted by supernova remnants (and other gamma transparent sources) . 
Namely, we calculate semi-analytically the (inverse) operator which relates an arbitrary gamma 
ray flux to the parent cosmic ray spectrum, without relying on any theoretical assumption about 
the shape of the cosmic ray and/or photon spectrum. We illustrate the use of this technique by 
applying it to the young SNR RX J1713. 7-3946 which has been observed by H.E.S.S. experiment 
during the last three years. Speciflc implementations of the method permit to use as an input either 
the parameterized VHE gamma ray flux or directly the raw data. The possibility to detect features 
in the cosmic rays spectrum and the error in the determination of the parent cosmic ray spectrum 
are also discussed. 

PACS numbers: 13.85.Tp, 96.50.sb, 98.70.Rz, 98.38.Mz 
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1 Introduction 

There is no doubt that the main part of cosmic rays (CR) till the knee is produced in the 
Milky Way pQ, and it seems fair to say that the favored site for CR production are the young 
supernova remnants (SNR). In fact, the turbulent gas of SNR is a large reservoir of kinetic 
energy [2] and this environment can support diffusive shock waves acceleration [3J. In recent 
times, great progresses have been made both in the observation and in the understanding of 
SNR. In particular, the new generation imaging Cherenkov telescopes, in particular H.E.S.S. [4J, 
allowed to observe the very high energy (VHE) gamma rays emitted by SNR, which are possibly 
generated by the decay of 7r° and t] produced by collision between the accelerated hadrons and 
the surrounding gas. It is not yet possible, however, to exclude that (part of) the observed 
radiation is produced by electromagnetic processes. In order to reach a definitive proof of the 
hadronic origin of the VHE gamma radiation, more detailed studies are needed. In this respect, 
new data at high (100 TeV or larger) and low {E^ ~ energies, improved theoretical 

modeling and possibly observations of VHE neutrinos (see e.g., [5]) will be extremely important. 

The hypothesis that VHE gamma radiation from young SNR originates from hadronic pro- 
cesses deserves the most serious attention and consideration. New and crucial observations are 
being collected and the hadronic origin seems to be favored for certain SNR, such as Vela Jr [B] 
and RX J1713. 7-3946 [TjJll In this paper we take the hadronic origin has a working hypothesis 
and we address the question of what we learn on SNR cosmic ray spectra from VHE 7— ray 
data. This question has a precise quantitative character and we answer it in the most direct 
way. Namely, we calculate semi-analytically the (inverse) operator which relates an arbitrary 
gamma ray flux to the parent cosmic ray spectrum, without relying on any theoretical assump- 
tion about the shape of the cosmic ray and/or photon spectrum. We then illustrate the possible 
applications of our method by considering the H.E.S.S. data of RX J1713. 7-3947 that reached 
an impressive accuracy in the energy range from 300 GeV to 300 TeV [9]. We remark that in 
this case (and, more in general, whenever the source shows non trivial spectral features) the 
approximation of power law distribution and the many techniques of calculations tailored to this 
assumption are not adequate. 

The plan of the paper is as follows. In Sect. [2] we formulate the problem and we obtain a 
general, analytical solution. In Sect. Owe consider possible applications of our results. First, we 
derive the parent cosmic ray flux of RX J1713. 7-3946 by using suitable parameterizations of the 
gamma ray data. Then, we extract the information directly from the observational data. This 

^ See also [S] for a recent analysis leading to different conclusions. 
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second technique requires fewer assumptions and allows to propagate the observational errors 
easily. However, when applied to noisy data it requires a sort of image processing (Gaussian 
smearing) to produce a reasonable result. In Sect. H] we summarize our results, putting emphasis 
on the possible applications of our method. 



2 How to invert the relation between the photon and the CR spectrum 
2.1 Formulation of the problem 

We assume that the VHE photon flux in SNR has a hadronic origin, i. e., gamma-ray are produced 
by a flux of high energy cosmic ray protons interacting with an hydrogen ambient cloud having 
density n. Inelastic proton-proton interactions result in the production of 7r° and //-mesons which 
subsequently decay producing gamma-rays. It is important to note that SNR are 'transparent 
targets' for cosmic rays, as can be understood by very simple estimates. The column density of 
the system is indeed much smaller than the TeV-proton and photon interaction lengths (Ap = 
rrtp/a ~ 40 gi/cw? and Xq ~ 60 gr/cm^), being: 

dz = n dl rrip = 1.5 x 10^^ gr/cm^ 
for n = 100 prot./cm^, dl = 3 pc. 

The possibly overestimated value for n corresponds to proton number density in a typical molec- 
ular cloud that could be associated to the SNR, and the distance dl is the one covered in 1, 000 yr 
at a speed of 3, 000 km/s. In other words, the proton and photon interaction probabilities (equal 
to dz/Xp and dz/X^ respectively) are 10~^ or smaller, so that proton multiple interactions and/or 
re-absorption of the produced photons are absent for the typical conditions of a young SNR. 

The gamma-ray flux produced on a detector placed at a distance R by cosmic ray 

protons interacting with a 'transparent' medium can be written a^: 



c r 



dEp dE^ 



where da^/dE^ is the inclusive cross-section for 7 production. Here n and drip/ dEp are the 
target hydrogen number density and the cosmic ray proton number density (per unit energy) 



^ This relation is valid if the CR momentum distribution is approximatively isotropic. If this assumption is 
removed one has to replace, here and in the following: 



1 dnp[r,Ep] dnp[r,Ep,n] 



47r dEp dQp dEp 



where dup/dilp dEp is the CR proton number density per unit energy and unit solid angle, n is the unit vector 
in the direction connecting the SNR to the detector and we have taken into account that the produced photons 
are almost collinear with CR protons. 
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respectively. Both of them depend on the position inside the SNR, indicated by the coordinate 
vector r. 

Next, we adopt the usual definition of the adimensional distribution function E'p], 
according to which: 



^7 p 

E' P 



where a is the total inelastic p-p cross-section given by [TU] : 

a[Ep] = 34.3 + 1.88 ln[Ep/lTeV] + 0.25 ln[Ep/lTeV]2 mb 



(2) 



(3) 



Hadronic interactions are affected by quite large uncertainties and independent calculations of 
[x,Ep\ may differ at the 20% level [12] . In this work, we use a simple analytic formula pre- 



sented in [To] (see appendix |A]) which describes the results of public available SIBYLL code [I3 
with a few per cent accuracy over a large region of the parameter space (x > 10~^ 
100 GeV). By using rel. ([2]), we can rewrite eq. ([T]) as: 



<I>^[E^ 



°° dEr, 



E^ 



pj 7 



^7 p 



(4) 



where we introduce the important quantity: 

c a[Ep] 



47ri?2 



J d^r n[r 



dnp[r,Ep] 
dE^ 



(5) 



The function $p[£^p] is the quantity which is constrained by and most directly related to the 
VHE gamma ray observations. It has the dimensions of a differential flux, and below we will 
use cm~2 s~^ TeV~^. In the following, we will refer to $p[i?p] as the effective cosmic ray flux 
from the SNR and we will show how this quantity can be calculated starting from the photon 
flux ^j[E^]. 

When comparing with theoretical predictions, one should note that the effective CR flux 
encodes not only the energy distribution of cosmic ray protons but also the (weak) energy 
dependence of the cosmic ray interaction probability in the SNR. We note, in fact, that eq. ([5]) 
can be rewritten as: 

cNa[Ep 



47ri?2 



Jp[Ep 



(6) 



where N - 
given by: 



/ d^r n[r] is the total amount of target hydrogen in the observed system and Jp[-E'p] 



is the weighted average the CR energy distribution in the SNR with a weight function propor- 
tional to the target hydrogen distribution. A part from the constant term cN/{AttR'^) (which 
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can be deduced if independent information on the SNR distance and on the amount of target 
hydrogen are given), the two functions ^'p[-Ep] and Jp[-Ep] differs by the energy-dependent factor 
l/cr[£'p]. It should be noted that the cross section is slowly varying with energy, so that the 
main spectral features of $p[-Ep] (such as the position and the sharpness of a cutoff/transition) 
always reflects the spectral features of Jp[£'p]. In the region where the spectrum is approxi- 
mated by a power law, the energy dependence of cr[Ep] accounts for a small difference between 
the spectral indices of ^p[-E^p] and $p[i?p] which can be easily quantified being of the order of 
din a/ din Ep ~ 0.07 in the energy range of interest. It is clear that the above formalism can 
be applied to any gamma transparent source (not only to SNR) where the VHE gamma are of 
hadronic origin, such as a SNR-molecular cloud association [TT] . 

The question of what can we learn on the effective cosmic ray flux $p from $^ boils down 
to the task of inverting an integral equation (eq. (jlj)). In the rest of Sect. [21 we argue that, 
assuming a quasi-scaling behavior of hadronic cross sections (accurate at the few percent level), 
this problem can be solved in good approximation by applying the differential operator: 



where a„ are appropriate numerical coefficients given in sect. 12.31 Stated more clearly, we claim 
that the approximate inverse of eq. (jl]), symbolically written as $^ = jF[$p], is simply given by 
$p ^ ^^[3'7]- This result does not rely on any theoretical assumption about the shape of the 
cosmic ray and/or photon spectrum. We thus provide a simple method to extract and study 
possible spectral features of the parent cosmic ray flux in the SNR directly from the observed 
VHE gamma radiation. 

Finally, to help the readers who are more interested in applications than in the formal 
derivation of our results, we anticipate that, to tackle the mathematical problem, it was necessary 
to introduce a number of definitions. In this paper, we indicate the natural logarithm of proton 
and photon energies with the symbols ^p and (see eq. (fTOl) in the next section). Moreover, it is 
convenient to multiply the cosmic ray and the photon fiuxes by a power laws in energy according 
to (/?p = $p • (£'p/lTeV)" and = $7 ■ (-S7/lTeV)", where a is an appropriate coefficient. We 
remark that the "fluxes" (fp and (p^ are related by a differential operator of the kind ([HD for any 
value of a. The values of the numerical coefficients a„ can be easily calculated for any adopted 
value for a (see eq. fl20|) and related discussion). 

2.2 Notation and quas/ scaZ/ng" approximation 

It is useful to perform some changes of variables and rewrite the integral (jl]) as: 
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—00 



(9) 
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where proton and photon energies are expressed in terms of the variables: 



Si = In 

The fluxes are rewritten in terms of: 



1 TeV 



2 = p,7 (10) 



V,[e,] = $p[e^p] e-p, ^^^^ 

and the integral kernel f[y,ep] is defined by: 

f[y,e,] = e[~y]-e'^y-F,[ey,e^^] (12) 

The Heaviside function is introduced in order to integrate over the entire real axis. The inclusion 
of the exponential factors in definitions flTT]) and f[T^ is, instead, motivated by the fact that the 
photon and CR proton spectra are expected to decrease approximatively as power laws in energy. 
For a proper choice of the parameter a, the functions ipi = ■ {Ei/1 TeV)° are, thus, expected 
to be nearly constant in the relevant energy range, highlighting the deviations from the pure 
power-law behavior. In the following, we find it convenient to set the value: 

a = 2.5 (13) 

which is particularly appropriate for the analysis of the H.E.S.S. gamma ray spectrum of the 
RX J1713. 7-3946 supernova remnant. We remark that the "fluxes" ipi[ei] and the integral kernel 
f[y, e] defined above are easily tractable, since they are square- integrable in Si and y respectively. 

In Fig. [H we show the behavior of the integral kernel f[y,ep] as a function of y for selected 
values of ^p. We see that the function f[y,ep] is peaked at y = —1.8 (which corresponds to 
Ep/E^ = exp[—y] ~ 6) and that it is marginally dependent on the assumed proton energy. In 
the following, we assume a quasi- s caling hehavioi for hadronic cross sections, i.e., we neglect the 
dependence of f[y,ep] on Sp and replace: 

f[y,ep]~^f[y]^f[y,el] (14) 

where is a fixed reference value for the proton energy. We have chosen the value = 6.9 
{i.e., E^ = 1000 TeV) which is appropriate to calculate the gamma ray flux in the energy region 
E'^ ~ 1 — 1000 TeV probed by the H.E.S.S. experiment. Our calculations and Fig. [T]show that 
the quasi-scaling approximation is adequate at the few percent level. 



2.3 Formal solution of the inverse problem 

In the quasi-scaling approximation, we can invert the relation between the effective CR proton 
flux and the photon flux (a Volterra integral equation of the first type) by a simple semi-analytical 
method which gives very precise results. We obtain, in fact, a convolution integral: 

/oo 
dsp ipp[ep] f[e^-ep] (15) 
-oo 
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Figure 1: The integral kernel f[y,£p] as a function of y = ln[E^/Ep] for selected values of the proton 
energy Ep = ln[£'p/lTeV]. We have chosen Sp = 0, 2.3, 4.6, 6.9 corresponding to Ep ~ 1, 10, 100, 1000 
TeV, respectively. 

which can be treated by using standard techniques, such as Fourier analysisjf] finding 

where f[k], (p^[k] and (Pp[k] are the Fourier transforms of the functions f[y], ^p-yle-y] and v^p[£^p] 
respectively. We note that the inclusion of the exponential factor in the definition of the <fi[ei] 
ensures that the Fourier transforms (p-ylk] and fp[k] exist0 The function: 

m . (17) 

defines, in Fourier space, the operator which inverts eq. f|T5l) . We see from Fig. [2] that Abs[/i[/c]] 
is fast increasing with \k\. This can be understood in simple terms by noting that the integral 
kernel f[y] has a half-width-half- maximum equal approximatively to 6y ~ 1.0 (see Fig. [1]). 
Correspondingly, the Fourier transform f[k] has a characteristic width 6k ~ l/(27r5?/) ~ 0.16 
and its inverse function h[k] has a sharp increase for \k\ > 6k. In physical terms, this has the 
important consequence that any feature in the photon spectrum on scales smaller than 6e^ < 6y 
will be greatly amplified in the parent CR proton spectrum (more discussion in the next section) . 

The behavior of h[k] at large k depends on the regularity of f[y] and its derivatives. In this 
case, we can expand h[k] to fifth order in a Taylor series: 

5 

h[k]c^Y.^J^^ (18) 
7=0 



^We calculate Fourier transforms according to the standard definition: (p[e] — J dk (p[k] exp[27riA:£]. 
* The functions (pi[ei] decrease exponentially for \ei \ oo provided that the differential energy spectra 
decrease slower that E^" at low energy and faster than E^°' at high energy. 
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Figure 2: Absolute value and argument of the function h[k] defined in eg. (17). 

where hj = (1/j!) d^h/dy\k=Q, with a few per cent accuracy in the relevant range \k\ < 2.5. We 
remark that the photon flux v^7[£-y] is sampled by H.E.S.S. experiment in bins 6e^ = 6E^E^ ~ 
0.2 or larger and, thus, only the region \k\ < l/(2fc^) = 2.5 carries physical information|j 

The expansion fllSp allows to express the parent proton spectrum as a function of the photon 
flux and its derivatives. By exploiting the properties of Fourier transforms one easily obtains: 

5 dim 
j=0 "^7 

where aj = hj/{2iTiy. The coefficients aj depend on the value of the parameter a adopted in 
eqs. ( ITT|T2l) . In our case (a = 2.5) the relevant coefficients are given by ao = 20.85, ai/ao = 
-2.336, aa/ao = 2.113, Og/ao = -0.9034, 04/00 = 0.1718 and a^/ao = -9.79 ■ lO'^. For a 
different choice a —>■ a — (3, the coefficients have to be replaced by: 

For instance, if we set P = a which corresponds to the particular situation considered in eq. 
[i.e., (fp = $p and ip-y = $^), we immediately obtain oq = .1148, ai = 2.390, 02 = 5.205, 
03 = 4.225, 04 = 1.031 and 05 = -0.2041. 

The above equations are the main results of this paper and, in the next section, we will 
discuss the possible apphcations to real data. Here, we note that rel. f|T9|) is remarkably simple 
if the photon spectrum can be approximated by a power law, i.e., $^ oc or, equivalently, 
if^ oc exp[(a— r)£:..y]. We obtain, in fact: v^p[£^p] = 3^[r] '/'^[e^p] with y[T] = Yfj=o (^j {a — Ty which 
shows that the ratio between the effective cosmic ray flux and the photon flux at a fixed energy 
is given by the function 3^[r] which only depends on the photon spectral index F. The function 
y[T] can be compared with the spectrum weighted moments Z[T] displayed in Fig. 5.5 of [H]. We 
obtain a good agreement by noting that, in the assumption of [H], one has Z[T] ~ r/(2 3^[r]). 



The physically significant range |fc| < 2.5 has been estimated by applying the sampling theorem. 
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Finally, we remind that the differential operator in the r.h.s. of eq. (fT9|) is the inverse of 
the integral operator in eq. (fT5|) . which is obtained in the quasi-scaling assumption for hadronic 
cross-section. One can go beyond this approximation by using eq. (fT9|) as zero-order solution 
and calculating perturbatively the effects of deviations from the quasi-scaling assumption. We 
used this approach to check that the corrections in the parent cosmic ray spectrum obtained 
from f|T9|) are small in comparison to the hadronic cross section uncertainties. We do not discuss 
the numerical implementation here to avoid unnecessary complications, but we provide the 
relevant details in the appendix [Bl 



3 Applications 

3.1 The young SNR RX J1713. 7-3946 

The RX J1713. 7-3946 Supernova Remnant has been observed by H.E.S.S. during three years 
from 2003 to 2005 [9]. The 7-ray spectrum obtained by combining the observation of the three 
years is shown in Fig. O The data extend over three decades, exploring the energy interval 
Er^ = 0.3 — 300 TeV. The energy resolution of the experiment is equal to about 20% and the 
photon spectrum is sampled in 25 bins 6e^ = 6E^/E^ ~ 0.2 plus three larger bins at high 
energyj^ 

The described data allow to obtain important conclusions, as discussed in [H]. First, they 
show that there is a significant emission at energy larger than 30 TeV, implying the existence 
of primary particles of at least that energy. Moreover, the data show a non trivial dependence 
on energy. In particular, there is a significant deviation from the simple power law behavior, as 
can be understood from Fig. [31 The solid line is the best fit power law spectrum with spectral 
index F = 2.32 that does not provides an acceptable fit of the data since x^/d.o.f. ~ 145.6/25 
{i.e., can be rejected at 9cr). 

From a theoretical point of view, one expects that the photon spectrum can be described by 
power law at low energy with a "cutoff' above a certain energy E^ related to to the properties 
of the primary particles acceleration mechanism. This kind of behavior is usually parameterized 
in the form of a broken power law (BPL): 



% = ^(^y^ U+fl^) ) [BPL case] (21) 



1/5^ -S(r2-ri) 

VlTeVy V" ' 

where Fi and F2 are the low and high energy spectral indices and S quantifies the sharpness of 
the transition from Fi and F2, or by an exponential cutoff (EC) with exponent f3: 



( -^7 A 



-r 



Ej 



[EC case] (22) 



^To help readability, in Figs. [3]l6] we use the logarithm to basis 10, denoted by log. 
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Figure 3: The j-ray spectrum of the RX J1713. 7-394-6 Supernova Remnant obtained by H.E.S.S. 
experiment. The black line represents the best fit to the data in the assumption of a power law behavior 
of the photon spectrum. The blue lines are obtained by assuming that the photon spectrum follows the 
BPL given by eq. 121\) in the assumption that the transition parameter is S = 0.6 (blue solid line) or 
S = 0.45, 0.75 (blue dotted lines). The red lines refers to the EC case, see eq. [2^). in the assumption 
that (3 = 0.5 (red solid line) or (3 = 1.0 (red dotted line). 

The value (3 ~ 0.5, which describes a relatively smooth cutoff, has been considered in [ini fT5] . 

The H.E.S.S. gamma ray spectrum of the RX J1713. 7-3946 supernova remnant has been fitted 
with a BPL with parameters Fi = 2.00 ± 0.05, F2 = 3.1 ± 0.2 and = 6.6 ± 2.2 obtaining a 
X^/d.o.f. ~ 29.8/23 (see Fig. [31 blue solid line). It should be noted that the sharpness parameter 
5" was kept fixed in the fit, with an adopted value equal to S* = 0.6. Different choices for S, 
however, are possible. As an example, equally good fits of the data are provided by the blue 
dashed line which corresponds to S* = 0.75, Fi = 1.97, F2 = 3.22 and Ec = 7.97 and by the blue 
dotted line which corresponds to 5* = 0.45, Fi = 2.03, F2 = 2.96 and Ec = 5.64. 

Alternatively, the H.E.S.S. data can be fitted with an EC with parameters F = 1.79 ± 0.06, 
Ec = 3.7 ± 1.0 and P = 0.5 as it is shown by the red solid line in Fig. [3l obtaining a x^/d.o.f. ~ 
34.3/24. For comparison, we also show with a red dotted line the best fit which obtained by 
assuming "pure" EC {i.e., /3 = 1.0). In this case, one obtains F = 2.04 ± 0.04, Ec = 17.9 ± 3.3 
and a slightly worse fit to the data x^/d.o.f. ~ 39.5/24 {i.e., a goodness of fit of 2.4%). 

As discussed in [3 [9] , the observed spectral shape seems to favor the hadronic origin. In the 
following, we assume the hadronic origin as a working hypothesis and we discuss what the data 
can tell us about the primary proton spectrum in RX J1713. 7-3946. 
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Figure 4: Left Panel: The effective CR flux $p[£^p] from the SNR RX J 171 3. 7-3946 obtained from 
BPL and EC parameterizations of the gamma-ray flux measured by the H.E.S.S. experiment. The blue 
lines are obtained from the best-fit BPL parameterizations with sharpness parameter S = 0.6 (blue solid 
line), S = 0.45 (blue dotted line) and S = 0.75 (blue dashed line). The red lines correspond to best-fit 
EC parameterization with (3 = 0.5 (red solid line) and (3 = 1.0 (red dotted line). Right Panel: The 
CR energy distribution J[Ep\ calculated from eq. ^ by assuming R = 1 kpc and N = 3.6 x 10^^. 

3.2 Using parameterized fluxes 

If we accept the BPL and EC parameterizations as reliable descriptions of the photon spectrum, 
we can calculate the effective CR flux from the SNR by simply applying eq. (fT9|) to the functional 
forms (12T]) and (122|) . The results of this procedure are shown in Fig. H] (left panel) where the 
blue lines are obtained from the BPL parameterizations of the photon flux, while the red lines 
refer to the EC case. We remind that the effective cosmic ray flux encodes not only the energy 
distribution of cosmic ray protons but also the (weak) energy dependence of the cosmic ray 
interaction probability in the SNR. For this reason we also show (right panel) the CR energy 
distribution in the SNR calculated according to eq. IQ. We assume R = 1 kpc and = 3.6 x 10^^ 
which corresponds to 300 solar masses of target hydrogen. This value is motivated if gamma ray 
emission is due to a molecular cloud-SNR association of the type proposed in [17] which seems 
to be consistent with the observations of NANTEN [18]; see [19] for a theoretical model. In this 
work, we do not aim to discuss the precise value of A^, which would only affect the normalization 
of the CR energy distribution. We focus instead on the CR spectral properties which are directly 
determined by the observed photon spectrum. We remark a few important points. 

Accuracy of the inversion. The obtained CR fluxes can be used in rel. f|T5l) in order to check 
the accuracy of the inversion method. In all cases, the re-calculated photon fluxes agree with the 
input photon flux {i.e., adopted in rel. f|T9|) ) at the level of few parts per thousand in the energy 
range E.y = 1 — 1000 TeV. This show that the differential operator on the r.h.s. of eq. (HM is the 
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inverse of the integral operator in eq. (fT5l) with very good accuracy, especially when compared 
with the uncertainties in the hadronic cross-section (at the ~ 20% level) or with the accuracy of 
the quasi scaling approximation (at the level of few percent or better). 

Cutoff /transition in the CR spectrum. The calculated CR spectra indicate, in all cases, that 
there is a significant number of protons at high energy. Protons should be efficiently accelerated 
up to an energy equal to about ~ 100 TeV, in order to explain the observed data. We see that 
the cutoff /transition region in the CR spectrum is in the energy region ii^p = 30 — 150 TeV. It 
is interesting to note that a photon flux with a smooth EC {(3 = 0.5) corresponds to an effective 
CR flux $p[ii^p] which is well described by the simple functional form exp[—E / E^] with 
r ~ 1.79 and cutoff energy E^ — 113 TeV, in reasonable agreement with the conclusion of |15j . 
The corresponding CR energy distribution Jp[-E'p] is well fitted by the same functional form, with 
the same cutoff energy and with F ~ 1.86. Compare with the discussion after eq. ([7j). 

Other features in CR and photon spectrum. We note that the differences between the CR 
spectra are much larger than the differences between the input photon fluxes. In the BPL case, 
the calculated CR spectra have a complex behavior in the energy range i?p = 3 — 30 TeV. 
Similarly, in the EC case the obtained curves differ substantially in the energy region Ep = 
30 — 150 TeV. This is not an artifact of the inversion method which is accurate at the level of 
few parts per thousand or better. It simply reflects the fact that any sharp feature in the photon 
flux is amplifled in the parent CR spectrum. In particular, the sharper is the transition/cutoff 
in the photon spectrum {i.e., the smaller is the S in eq. or the larger is the f3 in eq. (122]) ). 
the more complex is the behavior of the calculated CR flux. 

Dilution of spectral features. The previous point can be understood in terms of the properties 
of hadronic cross sections. It is basically related to the fact that the photon spectrum is not 
supposed to have any sharp feature if it is originated by hadronic processes. The integral kernel 
f[y] has, in fact, a characteristic width 6y ~ 1.0. Consequently, features in the CR spectrum 
on scales 6ep < 6y are washed-out by convolution f|T5l) . Conversely, if we observe features in 
the photon spectrum on scales 6e^ < 6y we are forced by rel. f|T9l) to postulate a complicated 
behavior of the parent CR proton flux which may be difficult or impossible to justify. As an 
example, BPL fits of the observational data with S < 0.4 correspond to parent CR spectra which 
become negative in the region Ep = 3 — 30 TeV and are, thus, not acceptable. 

A plausibility test for the hadronic origin assumption. The presence (or the absence) of features 
in the observed photon fiux on scales 6e^ < 6y may be used, in principle, as plausibility criterion 
to reject (or support) the hadronic origin of the observed 7-ray fiuxes|l| Interestingly, in the 

^ The energy resolution of the H.E.S.S. experiments {Se^ ~ 0.2) is sufficiently good, in principle, to test 
whether there is some sharp feature in the photon spectrum. 
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case of EC parameterization, the experimental data prefers a smooth cutoff (/3 ~ 0.5) which 
corresponds to a simple behavior of the primary photon flux. In the BPL case, our results 
shows instead that statistical errors are too large to arrive at any relevant conclusion about the 
sharpness of the transition. Thus, the fine structures of the primary CR spectra in the energy 
range ii^p = 3 — 30 TeV are not significantly constrained by the data. 



3.3 Using the raw data 

The differential operator on the r.h.s. of eq. (fT9l) is the inverse of the integral operator in eq. (fTSl) 
and it allows to obtain the CR spectrum directly from the photon flux, independently of any 
theoretical assumption. The CR spectrum, however, depends on high order derivatives of the 
photon flux which are generally known with bad accuracy and, moreover, introduce complicated 
correlations between the values of the CR flux extracted at two different energies. This makes 
difficult to infer the parent CR flux directly from noisy data and one could be tempted to 
conclude that a parameterization of the gamma ray flux is, in fact, necessary. In this section we 
propose a non-parametric procedure (based on Gaussian smearing) that avoids these difficulties 
and moreover permits to evaluate the error on the inferred CR spectrum. 



The 'smoothing' procedure. The relevance of the high order terms in rel. (fT9i) depends on the 
scale of the features that we probe in the CR spectrum. If we are interested in scales 6ep > 6, 
we can define the smoothed CR spectrum as follows: 



/oo 
de ipp[e] r[ep-e,6] 
-oo 



where: 



r[e,5] 



2n5 



exp 



e 



(23) 



(24) 



In Fourier space, this is equivalent to apply a Gaussian filter to fp[k] with a width equal to 
Ak = l/(27r5). We remind, for clarity, that (pp = <J>p [Ep/lTeV)'^'^ where the effective cosmic 
ray flux $p is defined in eq. ([5]). Equivalently, we can define the smoothed CR energy distribution 
by: 

/oo 
de jp[€]r[ep- e,S] (25) 
-00 



where Jp = x (i^p/lTeV)^'^ and Jp is given in eq. ([7]). By using eq. 
that: 

AttR^ 



it is possible to show 



(26) 



cNa[ep] 

with a few per cent accuracy, in the energy range of interest. We will then focus on the calculation 
of Ipp, showing that it can be simply estimated from observational data. 
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The 'smoothed' CR spectrum. Applying the differential operator of eq. (fT9|) . we find that the 
smoothed CR spectrum Tf^ is related to the photon flux by a convolution integral: 

/oo 
de^ ip^[e^] p[ep - 6^,6] (27) 
-oo 

where the convolving function p[e, 6] is given by: 

p[e,6]=r[e,6] ^A.e^ (28) 



2 = 



and the coefficients Ai are equal to: 



Ao-ao- — + 6—, --— + d— -15— , As-^-b— , 



(29) 



We can apply rel. fl27j) directly to the raw data, as it is explained in the following. We indicate 
with ifi ± Aipi the value of the photon flux measured in the i-th bin, centered at the photon 
energy Ei and covering the energy range (£i,inf, £^j,sup)- We approximate the photon flux by: 

ip4e^] = Y.ip,W,[e^] (30) 

i 

where VFi[^7] rectangular functions describing the various energy bins {i.e., Wi[ey] = 1 for 
^i.inf < < £i,sup and zcro elsewhere). We immediately obtain from eq. (!27|) the relation: 

^p[^P,^] = J2ViWi[ep,6] (31) 

i 

where: 



Wi 



[sp, 5]= de^ p [sp - e^, 5] ^ p [sp - e^, 5] Asi (32) 

and, in the last step, we have assumed that 6 3> A^j = ^j^sup — ^j,inf- Eq. fl?T]) gives the desired 
expression of the (smoothed) CR flux direcly from the gamma ray data. 

The smoothed CR spectrum is a linear combination of the observational values (pi of the 
photon flux. The functions Wj[£:p,(5] describe the contribution that each data point give to the 
reconstructed spectrum at a fixed energy ^p. The uncertainty in the CR spectrum can be easily 
evaluated by propagating linearly the observational errors Aipi obtaining: 



A^p[ep,(5] ^ ^J:^Avlw,[ep,SY ^^^^ 
'fp[£p,S] Y.i'fiWi[€p,5] 

Similarly, the correlation between the values of the CR flux at two different energies can be 
obtained by calculating: 

r Ek^vlwk[ep,5]wk[e'p,5] 
g\ep,e„,6\ = , = — , (34) 
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Figure 5: Left Panel: The (smoothed) CR spectrum from the RX J1713. 7-3946 SNR deduced from 
the raw data of the H.E.S.S. experiments. The solid line is obtained by continuing the j-ray spectrum 
at low and high energy with the best-fit BPL with sharpness parameter S = 0.6, while the dashed line is 
obtained by using the best-fit EC with [j = 0.5. The shaded area represents the observational uncertainty 
and is obtained by propagating H.E.S.S. observational errors. Right Panel: The (smoothed) CR 
energy distribution calculated from eg. h26\) with R = 1 kpc and N = 3.6 x 10^^. 

Application to the RX J1713. 7-3946 observations. Before applying the above relations to the 
H.E.S.S. data we have to choose the smoothing scale S. The choice of 6 is somewhat arbitrary 
and depends on the detector, on the quality of the observational data and on the problem under 
consideration. The H.E.S.S. detector has an energy resolution equal to 6e^ = 0.2 that suggests 
to adopt 6 ^ 0.2. Moreover, if we accept the hadronic origin assumption, we know that the 
hadronic interactions themselves introduce a scale, Sy ~ 1.0, below which the photon spectrum is 
not expected have features large enough to be significant with respect to observational errors. At 
the same time, we know that "noise" at these small scales is greatly amplified in the calculated 
CR spectrum. All this suggests to choose 6 = 6y = 1.0 and to focus our attention on the large 
scale features of the parent CR flux. 

Our final results are displayed in Fig. [51 In the left panel we show the smoothed CR spectrum 
which is obtained from the H.E.S.S. observational data of the RX J1713. 7-3946 SNR. In the 
right panel we show the smoothed CR energy distribution estimated according to eq. fl26l) with 
R = 1 kpc and = 3.6 x 10^^. We remind that the observational data cover an energy range 
equal to = 0.3 — 300 TeV. In this energy range we have described the photon spectrum 
according to eq. fl30|) . while at low < 0.3 TeV) and high energy {E^ > 300 TeV) we have 
continued the photon spectrum by using the best-fit BPL with sharpness parameter S = 0.6 (see 
previous section)!^ The shaded area describes the observational uncertainty in the smoothed CR 
spectrum which is obtained by propagating the errors in the observational data according to 

* Strictly speaking, one should know the photon flux at all energies in order use the rel. (j3ip . 
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Figure 6: Correlation between the values of the (smoothed) proton spectrum obtained at two different 
energies. Light colors correspond to values of the correlation index close to 1, while dark colors corre- 
spond to values close to -1. The displayed contours corresponds to correlation index equal to 0.5 and to 
-0.5. 

eq. (l33l) . One sees that the error is less than 10% at low energy and remains smaller than 20% 
for i?p < 300 TeV. The correlation between the values of the CR flux at two different energies, 
evaluated as in eq. is shown in Fig. [61 

In order to estimate the systematic uncertainty introduced by the ignorance of the high 
and/or low energy behaviour of the photon flux, we also show with a dashed line the smoothed 
CR spectrum which is obtained by continuing the photon spectrum at high and low energy with 
the best-fit EC with (3 = 0.5. The difference between the solid and the dashed line is smaller than 
the observational error for Ep = 3 — 1000 TeV, showing that the proton spectrum, in this energy 
range, is directly constrained by observational data. For Ep < 3 TeV, instead, the two lines 
behave differently indicating that the systematic uncertainty due to extrapolation is relevant. 
In this respect, the small bend of the effective CR flux in the energy region 1 — 10 TeV, at the 
moment, does not seem to be fully signicative. The existence of such a bend would amount to 
an important physical information on the acceleration mechanism (see e.g., [IS]). Thus, it will 
be interesting to collect new dcitci cit energy lower than E^ = 0.3 TeV to assess its significance. 

In conclusion, the displayed results show that the large scale features of the CR spectrum in 
the energy range Ep = 3 — 300 TeV are well constrained by the observational data. The effective 
CR flux is roughly described by power-law with spectral index F = 1.7 — 2 at low energy with a 
cutoff /transition region between ii^p = 30 — 100 TeV. This conclusion is also consistent with the 
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one obtained in Sect. 13.21 using parameterized photon fluxes. 



4 Summary 

In this work we assumed the hadronic origin of the gamma radiation emitted by SNR and we 
addressed the question of what can be learned on the SNR cosmic ray spectrum from VHE 7-ray 
data. We summarize here our conclusions: 

i) The main result is contained in eq. fll9l) . This equation shows that, in the approximation 
of quasi-scaling behavior of the hadronic cross sections]^ the effective CR spectrum defined in 
eq. can be obtained by applying a simple differential operator to the observed photon flux. 
This results does not rely on any theoretical assumption about the shape of the proton and/or 
photon spectrum. It can thus be applied to sources which show non trivial spectral features such 
SNR RX J1713. 7-3946 for which, instead, the commonly adopted approximation of power law 
distribution and the many techniques of calculations tailored to this case are not adequate. 

ii) We have emphasized that the presence (or the absence) of sharp features in the photon 
spectrum can be used as plausibility criterion to reject (or to support) the assumption that the 
observed radiation has a hadronic origin. The basic point is that the hadronic processes (con- 
volved with the parent CR flux) lead to a characteristic energy scale below which the produced 
photon flux is expected to be featureless (see discussion in Sect. 13. 2p . 

Hi) Specific implementations of our method permit to calculate the parent CR spectrum either 
from parameterized VHE fluxes or directly from raw data (see eq. flSTl) ). This second approach 
requires fewer theoretical assumptions and allows to propagate the observational errors easily 
(see eqs. fl33ll34p ). However, when applied to noisy data, it requires a sort of image processing 
(Gaussian smearing, see eq. (123|) ) to produce reasonable results. 

iv) We have applied our method to the young SNR RX J1713. 7-3946 which has been observed 
by the H.E.S.S. experiment during the last three years. We have calculated the CR spectrum both 
from parameterized photon fluxes and directly from the raw data. The results are summarized in 
Fig. m and Fig. O These figures demonstrate that the observational data constrain well the main 
features of CR flux in the energy range E'p ~ 3 — 300 TeV; they give, instead, a poor information 
outside this range, and cannot significantly test the fine structures of the CR spectrum. As a 
final result, we conclude that the effective CR flux from SNR RX J1713. 7-3946 is well described 
by power-law with spectral index F = 1.7 — 2 at low energy with a cutoff/transition region 
between = 30 - 100 TeV. 

^The "quasi scaling" assumption defined in eq. ITU) is accurate at the few per cent level, see Fig. [1] and can 
be improved as discussed in the appendix [Bl 
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A The functions F^[E^/E^^E^ and f[y 

In order to give a self-contained discussion, we give here the function F^[E^ / E^, E^] obtained 
in [To] and used in this paper. We have: 

\l + k^x^^{l-xP-<)) \]^x\ ~ l-xl^-'" 1 + k^x^'<{l - xP'<) J 

where x = E^/Ep. The parameters B^, (3^, and k^ depend only on the energy of proton and are 
given by: 

B^ = 1.30 + 0.14£p + 0.011 4, 

f3, ^ ^ 



1.79 + 0.11 ep + 0.008 el 



0.801 + 0.049 Ep + O.OUel 



where = ln[i?p/l TeV]. The function f[y] is obtained by applying eqs. (|T2|) and (IT^ and is 
simply given by: 

/[!/] = ^[-y]■e"^•F,[e^lPeV] 



B Improving the quasi scaling approximation 

Let us begin by writing the integral equation eq. ([5]) in abstract terms: 

$,[E,]=^[$p][i?,] (35) 

We showed that a solution of the integral equation $^ = J^o['^'p] in the quasi-scaling approxi- 
mation F[x,Ep] — > F[x,Ep], accurate at the level of few parts per thousand or better, is given 
by 

$p = V[%] (36) 
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In operator terms, we write DJ^q = !FqD = 1. For a given gamma-ray flux $^ we can evaluate 
the goodness of a certain approximation of the proton 'flux' <I>p from the difference between the 
assumed gamma-ray flux and the gamma-ray flux re-calculated from the proton 'flux': 

i - (37) 

In our case, we evaluate the goodness of the quasi-scaling approximation by using eq. (!36|) for 
$p. Assuming that $^ is a broken power-law, we flnd that ^ is smaller than 10% in the range 
of energies from 100 GeV to 1 PeV, that is already sufficient for our purposes. It is however 
possible to improve the approximation for the proton flux as follows: 

$p = P[$,(l-0] (38) 

where ^ is calculated using eq. fl36l) l^ Repeating the procedure (namely: assuming again the 
broken power-law distribution for $^ and plugging eq. fl38l) into eq. fl37|) ) in order to test the 
approximation, the newly calculated ^ is smaller than 0.5% in the range from 100 GeV to 1 PeV. 



i"The formal derivation is simple: From <I> = = {To + SJ^)[^] we get To[i'] = <f> - Applying V we 

find \E' = — (5jr[\I>]], that can be improved iteratively. If in the r.h.s. of the last equation we use 4' = I?[<1>] 
{i.e., eq. ([55)) ) we get eq. (1551) simply applying the definition of eq. ([57)) . 
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